import numpy as np
import matplotlib.pyplot as plt

Fe = 4  #fréquence d'échantillonnage (Hz)
Te = 1 / Fe  #période d'échantillonnage (s)
duration = 10  #durée (s)

#vecteur temporelle
t = np.arange(0, duration, Te)

#definition du signal s(t)
s_t = 3 + np.sin(2 * np.pi * t)

plt.figure(figsize=(10, 4))
plt.plot(t, s_t, label="s(t) = 3 + sin(2πt)")
plt.xlabel("Temps (s)")
plt.ylabel("Amplitude")
plt.title("Signal s(t) sur 10 secondes")
plt.legend()
plt.grid()
plt.show()

from scipy.signal import lfilter

help(lfilter)

Help on function lfilter in module scipy.signal._signaltools:

lfilter(b, a, x, axis=-1, zi=None)
    Filter data along one-dimension with an IIR or FIR filter.

    Filter a data sequence, `x`, using a digital filter.  This works for many
    fundamental data types (including Object type).  The filter is a direct
    form II transposed implementation of the standard difference equation
    (see Notes).

    The function `sosfilt` (and filter design using ``output='sos'``) should be
    preferred over `lfilter` for most filtering tasks, as second-order sections
    have fewer numerical problems.

    Parameters
    ----------
    b : array_like
        The numerator coefficient vector in a 1-D sequence.
    a : array_like
        The denominator coefficient vector in a 1-D sequence.  If ``a[0]``
        is not 1, then both `a` and `b` are normalized by ``a[0]``.
    x : array_like
        An N-dimensional input array.
    axis : int, optional
        The axis of the input data array along which to apply the
        linear filter. The filter is applied to each subarray along
        this axis.  Default is -1.
    zi : array_like, optional
        Initial conditions for the filter delays.  It is a vector
        (or array of vectors for an N-dimensional input) of length
        ``max(len(a), len(b)) - 1``.  If `zi` is None or is not given then
        initial rest is assumed.  See `lfiltic` for more information.

    Returns
    -------
    y : array
        The output of the digital filter.
    zf : array, optional
        If `zi` is None, this is not returned, otherwise, `zf` holds the
        final filter delay values.

    See Also
    --------
    lfiltic : Construct initial conditions for `lfilter`.
    lfilter_zi : Compute initial state (steady state of step response) for
                 `lfilter`.
    filtfilt : A forward-backward filter, to obtain a filter with zero phase.
    savgol_filter : A Savitzky-Golay filter.
    sosfilt: Filter data using cascaded second-order sections.
    sosfiltfilt: A forward-backward filter using second-order sections.

    Notes
    -----
    The filter function is implemented as a direct II transposed structure.
    This means that the filter implements::

       a[0]*y[n] = b[0]*x[n] + b[1]*x[n-1] + ... + b[M]*x[n-M]
                             - a[1]*y[n-1] - ... - a[N]*y[n-N]

    where `M` is the degree of the numerator, `N` is the degree of the
    denominator, and `n` is the sample number.  It is implemented using
    the following difference equations (assuming M = N)::

         a[0]*y[n] = b[0] * x[n]               + d[0][n-1]
           d[0][n] = b[1] * x[n] - a[1] * y[n] + d[1][n-1]
           d[1][n] = b[2] * x[n] - a[2] * y[n] + d[2][n-1]
         ...
         d[N-2][n] = b[N-1]*x[n] - a[N-1]*y[n] + d[N-1][n-1]
         d[N-1][n] = b[N] * x[n] - a[N] * y[n]

    where `d` are the state variables.

    The rational transfer function describing this filter in the
    z-transform domain is::

                             -1              -M
                 b[0] + b[1]z  + ... + b[M] z
         Y(z) = -------------------------------- X(z)
                             -1              -N
                 a[0] + a[1]z  + ... + a[N] z

    Examples
    --------
    Generate a noisy signal to be filtered:

    >>> import numpy as np
    >>> from scipy import signal
    >>> import matplotlib.pyplot as plt
    >>> rng = np.random.default_rng()
    >>> t = np.linspace(-1, 1, 201)
    >>> x = (np.sin(2*np.pi*0.75*t*(1-t) + 2.1) +
    ...      0.1*np.sin(2*np.pi*1.25*t + 1) +
    ...      0.18*np.cos(2*np.pi*3.85*t))
    >>> xn = x + rng.standard_normal(len(t)) * 0.08

    Create an order 3 lowpass butterworth filter:

    >>> b, a = signal.butter(3, 0.05)

    Apply the filter to xn.  Use lfilter_zi to choose the initial condition of
    the filter:

    >>> zi = signal.lfilter_zi(b, a)
    >>> z, _ = signal.lfilter(b, a, xn, zi=zi*xn[0])

    Apply the filter again, to have a result filtered at an order the same as
    filtfilt:

    >>> z2, _ = signal.lfilter(b, a, z, zi=zi*z[0])

    Use filtfilt to apply the filter:

    >>> y = signal.filtfilt(b, a, xn)

    Plot the original signal and the various filtered versions:

    >>> plt.figure
    >>> plt.plot(t, xn, 'b', alpha=0.75)
    >>> plt.plot(t, z, 'r--', t, z2, 'r', t, y, 'k')
    >>> plt.legend(('noisy signal', 'lfilter, once', 'lfilter, twice',
    ...             'filtfilt'), loc='best')
    >>> plt.grid(True)
    >>> plt.show()

#coeffs du filtre moyenneur (fenêtre de 4 points)
coeffs_moyenneur = np.ones(4) / 4

#application du filtre moyenneur sur signal s(t)
s_moyenne = lfilter(coeffs_moyenneur, 1, s_t)

plt.figure(figsize=(10, 4))
plt.plot(t, s_t, label="Signal original s(t)")
plt.plot(t, s_moyenne, label="Signal filtré (moyenneur 4 points)", linestyle='--', color='r')
plt.xlabel("Temps (s)")
plt.ylabel("Amplitude")
plt.title("Filtre moyenneur 4 points")
plt.legend()
plt.grid()
plt.show()

#def du filtre dérivateur : différencier le signal par différence finie
coeffs_derivateur = [1, -1]

#application du filtre dérivateur au signal s(t)
s_deriv = lfilter(coeffs_derivateur, 1, s_t) / Te  #divise par Te pour obtenir l'approximation de la dérivée

plt.figure(figsize=(10, 4))
plt.plot(t, s_t, label="Signal original s(t)")
plt.plot(t, s_deriv, label="Signal dérivé (filtre dérivateur)", linestyle='--')
plt.xlabel("Temps (s)")
plt.ylabel("Amplitude")
plt.title("Filtre dérivateur")
plt.legend()
plt.grid()
plt.show()

#calcul de la dérivée théorique
#dérive théorique de 3*sin(2pi*t) est s'(t) = 2pi*cos(2pi*t)
s_deriv_theorique = 2 * np.pi * np.cos(2 * np.pi * t)

plt.figure(figsize=(10, 4))
plt.plot(t, s_t, label="Signal original s(t)")
plt.plot(t, s_deriv, label="Dérivée par filtre dérivateur", linestyle='--')
plt.plot(t, s_deriv_theorique, label="Dérivée théorique s'(t) = 2πcos(2πt)")
plt.xlabel("Temps (s)")
plt.ylabel("Amplitude")
plt.title("Comparaison entre la dérivée théorique et la dérivée par filtre")
plt.legend()
plt.grid()
plt.show()

from scipy.io import wavfile

import plotly.graph_objects as go

#load file
sample_rate, signal = wavfile.read("Note_Mi.wav")
#normaliser si le signal est en entier
signal = signal / np.max(np.abs(signal))

plt.figure(figsize=(10, 4))
plt.plot(np.arange(len(signal)) / sample_rate, signal)
plt.xlabel("Temps (s)")
plt.ylabel("Amplitude")
plt.title("Signal de la note Mi")
plt.grid()
plt.show()

/var/folders/2l/zssv4k554m3gbdphlp5n6dtm0000gn/T/ipykernel_1073/1160791914.py:2: WavFileWarning: Chunk (non-data) not understood, skipping it.
  sample_rate, signal = wavfile.read("Note_Mi.wav")

fig = go.Figure()

#ajout du signal audio à la figure
fig.add_trace(go.Scatter(x=np.arange(len(signal)) / sample_rate, y=signal, mode='lines', name='Signal Audio'))

#config graph
fig.update_layout(
    title='Signal Audio Interactif',
    xaxis_title='Temps (s)',
    yaxis_title='Amplitude',
    hovermode='x',  #affiche valeurs en fonction de la position de la souris
)

fig.show()

from scipy.fft import fft

#calcul de la TFD et module
signal_fft = fft(signal)
N = len(signal)
frequencies = np.fft.fftfreq(N, 1/sample_rate)
magnitude = np.abs(signal_fft) / N  #module de la TFD normalisé
phase = np.angle(signal_fft) #argument de la TFD normalis

print(signal_fft)

[-12.73790814-0.j          -5.58583009+2.49197276j
   0.88471877-2.93342473j ...  -9.95848843+0.5295149j
   0.88471877+2.93342473j  -5.58583009-2.49197276j]

#spectre monolatéral
positive_freqs = frequencies[:N//2]
positive_magnitude = magnitude[:N//2]
positive_phase = phase[:N // 2]

plt.figure(figsize=(12, 6))

#spectre de magnitude
plt.subplot(2, 1, 1)
plt.plot(positive_freqs, positive_magnitude)
plt.xlabel("Fréquence (Hz)")
plt.ylabel("Amplitude")
plt.title("Spectre de magnitude (monolatéral) du signal")
plt.grid()

#spectre de phase
plt.subplot(2, 1, 2)
plt.plot(positive_freqs, positive_phase, color='orange')
plt.xlabel("Fréquence (Hz)")
plt.ylabel("Phase (radians)")
plt.title("Spectre de phase (monolatéral) du signal")
plt.grid()

plt.tight_layout()
plt.show()

#ajouter bruit blanc
noise_amplitude = 0.2  #ajustez selon le niveau de bruit désiré
noisy_signal = signal + noise_amplitude * np.random.normal(size=len(signal))

plt.figure(figsize=(10, 4))
plt.plot(np.arange(len(noisy_signal)) / sample_rate, noisy_signal)
plt.xlabel("Temps (s)")
plt.ylabel("Amplitude")
plt.title("Signal bruité (note Mi avec bruit blanc)")
plt.grid()
plt.show()

from scipy.signal import butter, lfilter

#paramètres du filtre Butterworth (ordre et fréquence de coupure)
ordre_iir = 5  #ordre du filtre
fc = 150  #fréquence de coupure en Hz (ajustez en fonction du spectre)

#conception du filtre passe-bas
b, a = butter(ordre_iir, fc / (sample_rate / 2), btype='low')

#filtrage du signal bruité
filtered_signal_iir = lfilter(b, a, noisy_signal)

plt.figure(figsize=(10, 4))
plt.plot(np.arange(len(filtered_signal_iir)) / sample_rate, filtered_signal_iir)
plt.xlabel("Temps (s)")
plt.ylabel("Amplitude")
plt.title("Signal filtré avec un filtre IIR")
plt.grid()
plt.show()

from scipy.signal import firwin, convolve

#paramètres du filtre FIR
ordre_fir = 50  #ajuster pour un résultat similaire au filtre IIR
fc_fir = fc  #utiliser la même fréquence de coupure que pour l'IIR

#conception du filtre FIR passe-bas
fir_coeffs = firwin(ordre_fir + 1, fc_fir / (sample_rate / 2))

#application du filtre FIR
filtered_signal_fir = convolve(noisy_signal, fir_coeffs, mode='same')

plt.figure(figsize=(10, 4))
plt.plot(np.arange(len(filtered_signal_fir)) / sample_rate, filtered_signal_fir)
plt.xlabel("Temps (s)")
plt.ylabel("Amplitude")
plt.title("Signal filtré avec un filtre FIR")
plt.grid()
plt.show()

TD6 Traitement de Signal¶

EXERCICE I¶

Q1.1¶

Q1.2¶

Q1.3¶

Q1.4¶

Q1.5¶

Conclusion¶

EXERCICE II¶

Q2.1¶

Q2.2¶

Q2.3¶

Q2.4¶

Q2.5¶

Q2.6¶

Q2.7¶

Q2.8¶

Q2.9¶